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HIGH ORDER EXPLICIT SYMPLECTIC INTEGRATORS FOR 
THE DISCRETE NON LINEAR SCHRODINGER EQUATION. 

JEHAN BOREUX, TIMOTEO CARLETTI, AND CHARLES HUBAUX 



Abstract. We propose a family of reliable symplectic integrators adapted to 
the Discrete Non— Linear Schrodinger equation; based on an idea of Yoshida 1121 
we can construct high order numerical schemes, that result to be explicit meth- 
ods and thus very fast. The performances of the integrators are discussed, 
studied as functions of the integration time step and compared with some non 
symplectic methods. 



under periodic boundary conditions, {p{x + L, t), q{x + L, t)) = {p{x, t), q{x, t)) for 
all X and t, is a widely studied multiparticles system subjected to nonlinear inter- 
actions, that can be used to model several relevant physical phenomena [9J, ranging 
from optics to solid state and atomic physics, for instance Bose-Einstein conden- 
sates. 

Often scientists, once modeling physical nonlinear phenomena, impose the con- 
dition = iq{t) for all t, in this case the system reduces to the standard cubic 
NLS equation. We hereby adopt the viewpoint of [5] where the conjugacy relation 
between p{t) and q{t) is not imposed a priori. We will nevertheless show that if 
this relation holds at t = then it will be preserved by our numerical integration 
scheme. 

Most properties of NLS are related to the asymptotic behavior of its solutions, 
there is thus a strong need for integration schemes, allowing large step sizes, to eas- 
ily cover large time spans, without degrading the efficiency of the numerical method 
hence avoiding the introduction of spurious outcome in the numerical simulations. 
These goals can be achieved using symplectic integrators, that are specifically de- 
signed to preserve the energy and possibly other first integral of the system. Our 
analysis will be performed on the following one dimensional discretization of the 
Eq. ([1]), DNLS for short [5]: 



1. Introduction 



(1) 




(2) 




Where a + sign represents the defocusing case hereby considered, a — sign represents the 
focusing one and where a denotes the complex conjugate of the complex number a. 
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using the periodic boundary conditions qi+N = qi and pi+N = Pi for all / G 
{1,...,A^}. Let us observe that the solutions of the DNLS arc 0{h'^) close to 
the solutions of the PDE ([T|), thus in the limit h 0, the former converge to the 
solutions of the latter. 

The integration method hereby proposed is based on the development of an 
idea introduced by Yoshida ^ 12i . The strong improvement allowed by this method 
with respect to other ones available in the literature, relies on the fact that we 
are able to provide symplectic integrators, that are explicit ones, thus very fast, 
and whose order can be easily as high as the eighth one. Other methods based on 
the construction of generating functions, see for instance [TJ [7] are not particularly 
suitable because the Hamiltonian function describing the DNLS cannot be split in 
the form H(p,q) — T{p) + V{q), namely it is not a potential Hamiltonian system. 
This implies in fact that the obtained numerical integration schemes are implicit 
methods, and thus a large amount of CPU time is used to compute the new position, 
or the new momentum, using some Newton-like method. Let us observe that Gauss' 
methods [4 , namely symplectic versions of Runge-Kutta, are also implicit ones and 
thus they suffer of the above mentioned limitations. 

Finally let us observe that, because the system cannot be put in the form of a 
perturbation of an integrable system, the method proposed in 8 is no more useful 
here: one can easily get a second order method and with some additional work also 
a fourth order method TT, but no higher orders are obtainable. 

The paper is organized as follows. In the next section we briefly recall the 
Hamiltonian structure of the non linear Schrodinger equation on a one dimensional 
lattice; then Section[3]will be devoted to the presentation and to the construction of 
the symplectic integrator, whose properties will be numerically studied in Sectional 
Finally we sum up and draw our conclusions in Section [5l 



The system ([T]) will be studied imposing the standard spatial discretization, see 
for instance [5], thereby called Diagonal NLS, that reads: 



Observe that h plays the role of spatial discretization parameter and thus the 
first terms on the right hand sides stem from the discretized second order spatial 
derivatives. 

Let us remark that one could be interested in studying directly the system ([3]) as 
a relevant model of non-linear interaction on a discrete one dimensional of lattice. 

One can easily prove that the DNLS can be cast into the Hamiltonian formalism 
using the following Hamilton function 



1=1 '- 

and moreover the variables (p;, g;) G C" do satisfy the standard Poisson equations 
{PhPm} ^ {qi,qm} and {p/,g„i} = (5i,,„ V/, m e {1, . . . , iV} . 



2. The non linear Schr5dinger equation on a lattice 



(3) 




using the periodic boundary conditions: 

(4) qi^^=qi and pi+N ^ Pi \fl E {1, . . . , N} . 
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Remark 2.1. Let us observe that the DNLS system ([3]) possesses another first 
integral independent of the energy, namely: 

N 

(6) I{p,q) ^^piqi , 

1=1 

that under the assumption p(t) ~ q{t), is usually called the mass of the system. 

The aim of this work is to define a family of symplectic integrators based on 
the Yoshida symplectic scheme [12] adapted to the DNLS system ^ and to study 
their numerical properties in terms of preserved quantities, CPU time needed and 
accuracy of the results. 

Let us observe that system ([3]) is not quasi-integrable, i.e. it cannot be decom- 
posed as the sum of an integrable one and a "small perturbation" , thus we cannot 
use the high order symplectic schemes SABA^ or SBABn proposed by Laskar and 
Robutel [5], whose accuracy strongly relies on the smallness of the perturbation. 

3. A SYMPLECTIC SCHEME 

For a sake of completeness, let us briefly recall the integration method proposed 
in [12j to numerically reconstruct the orbit of an Hamiltonian systems H{p,q). 

Given any initial condition (p(0),g(0)) and an integration time span [0,r], we 
proceed by decomposing the integration interval into small pieces of fixed size t. 
The method results thus in a fixed step size integration method. Then, in each 
small interval, we approximate the time r flow of the Hamiltonian systems 

=-% 

by a composition of basic symplectic flows of the form exp (tcjLj^) and exp (TdjLg), 
where A{p,q) and B{p,q) are two suitable functions providing a decomposition of 
the Hamilton function, i.e. H = A + B, {cj,dj) suitable constants to achieve the 
wanted order of precision of the integration scheme and exp (tLh) is a shorthand 
notation to denote the flow at time t of the Hamilton function H. More precisely 
we are looking for 
(7) 

exp (tLa+b) = exp (rciLyi)oexp {TdiLB)o- ■ -oexp (TCfcLA)oexp {TdkLB) + 0{T") , 

for some positive integers k and n. 

The method is particularly efficient once the maps exp (rcjLA) and exp (TdjLs) 
are explicitely computable, which requires A and B to be simple enough. In par- 
ticular this is true whenever A and B depend only upon one group of canonical 
variables and thus the Hamiltonian is of the so-called potential form. Of course 
this is a strong requirement that cannot always be achieved by suitable change 
of coordinates. In the following we propose the decomposition of the Hamilton 
function ([5]) given by 

(8) Aip^q)^^±phf and = f i^^±l^M(^i±l^ . 

1=1 1=1 
Let us observe that even if A and B both depend on all the canonical variables, we 
are able to explicitely compute exp(TLyi) and exp(TLB) (see ? 13.11 and 13. 2p and 
thus to propose a completely explicit method. 
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3.1. Computation of exp (tLa)- The equations of motion of the system with 
"Hamilton function" are given by: 



(9) 



\/le{i,...,N} 



\pi ^ -2ipfqi 



from which it trivially follows that piqi = C; is a first integral for all I, in fact: 

rCiiflow m = piqi + piqi 



Hence Q simplifies into 



(10) 



yie{i,...,N} 



-2ipfqf + 2ipfqf - . 
\pi = -2ipiCi 



ji ='2iqiCi, 

whose solution with initial datum (p;(0), g;(0));=i,...,Ar is for all t: 



(11) 



yie{i,...,N} 



\pl(t) 

[qi{t) =e 

Finally the time t flow of A is given by: 

(12) = (e-2»fl«l>l,...,e^2^PiV9ivr^^^g 

where ^ denotes the transposed of the vector 



-2iCit 



2iCit 



PliO) 
9/(0). 



,qN) 

'qi,...,e- 



3.2. Computation of exp [tLb)- The equations of motion corresponding to the 
B part of the Hamiltonian function are: 



(13) 



We{i,...,iV} 



[qi = w: (qi+i ~ '^qi + qi-i) , 

that can be cast in a compact form using the periodic boundary conditions (|4]) by 
introducing the circulant matrix (2j M obtained by cyclically permute the TV-vector 
(-2,1,0,. ..,0,1): 



(-2 1 
1 -2 1 
1 -2 



(14) 



M 







1\ 










1 -V 



\l ... . 

that is the linear diagonal system 

fp- =^i,Mp 
\il ^^Mq. 

Thus formally the time r flow of the Hamilton function B is given by: 



(15) 



(16) {P,q'f = e^^Hp,q)'^ = {e-^^''p,el^^'''q)^. 

The circulant matrices have some useful properties that allow one to explicitely 
compute their eigenvalues, eigenvectors and hence their exponential. 
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Proposition 3.1 (Circulant matrix). Let M he the circulant matrix given by p4|) 
and let p = e^^^*/^, then 

(1) The eigenvalues of the matrix M are given by: fij = —2 + + p^^^^^\ 
j^Q,...,N~l; 

(2) For j — 0, . . . , N ~ 1, the eigenvector associated to pj is given by 



(3) Let W be the matrix whose columns are the eigenvectors Wq, . . . ,wpf-i, 
then W is unitary, namely WW''' — W'^W = 1^?. Where ^ denotes the 
transposed complex conjugated matrix; 

(4) The exponential of M is given by: 

N 



e^' = We''''"-W^ ^ m e {1, . . . , iV} (e*^/™ = T7 E 



Uk-i (i-m)(fc-l) 
k=l 



being Mdiag the diagonal matrix with the eigenvalues on the diagonal. 

Proof. Point (1) can be proved by a direct computation, for all j e {0, . . . , — 1}, 
one has: 



1 / \ 1 

\ : ] 



Vl+pj("-2)_2p^(' 



pJ(p-^-2+pJ) 



observing that p^ = 1, we thus have 



_2 + p^+p^/ I' 

' = 7^ I : I = 



Let us point out that the complex vectors (u'j)j=o,...,Af-i form an orthonormal 
set for (with the standard complex scalar product): 

N 1 ^ 1 ^ 

m— 1 m— 1 m— 1 

Let US also remark that the eigenvalues pj are actually real: 

Pj = -2 + pi + /T^'f^-i) =-2 + p> + p-^ ^-2 + 2 cos(27rj/iV) , 

and moreover the eigenvalues coincide in pairs, more precisely there are {N — l)/2 
distinct eigenvalues, if N is odd, and [N — 2)/2 ones if N is even: pj — pk if (j + k) 
mod iV = 0. 

To prove that W is unitary, let us denote by X = WW^ and compute the element 
J, k of such matrix. By definition 

AT N N ^ N 

1=1 1=1 1=1 1=1 

thus X = Ipf. A similar computation can be done for W''W. 

Finally let us observe that MW = WM^iag, in fact one trivially has 

MW = {poWol ■ . ■ \pN-lWN-l) = Mdiag, 
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and because W is unitary, W ^ = T4^^, thus M = WMdi^gW^ . This imphes that 
the exponential can be computed as fohows: 

The element /, m of e*^ is thus given by: 

N N ^ N 



1 ^ 

fc=i 

□ 



One can thus explicitely compute the flow of B given by (fH 
Corollary 3.2. T/ie time t flow of the Hamilton function B can be rewritten as: 

jp' 



(17) 



or explicitely for all I {1, . . . , N} 



' = We'i^'^^''""'W'^c 



(18) 



Pi — N l^m=\ l^k=\ e- ^ P Pn 



Let us observe that (jlSp can be rewritten in compact vector form in a way inspired 
by the solution of linear ODE, namely as a linear combination of eigenvectors, as 
follows: 



(19) 



= e "■'"''< Wo, p(0) > Wo H he ^^^""^ < WAr_i,p(0) > wtv-i 

= 6^^^° < wo, 9(0) > wo + • • • + e""?^''''"-^ <WN~i,q{0) >WN-i. 



3.3. The integrator. Once we have the explicit maps exp(TL^) and exp(rLB) 
one can construct the basic second order symplectic scheme Stormer-Verlet/Leap 
Frog: Y2{t) — ei^^e'^^^e^'^-*. Then Yoshida proved [T^ that one can find explicit 
suitable coefficients xi — 2-2^/3 ^0 = — 2^/'^a:i, such that the composition 

(20) Yi{T) = Y2 (xiT) o Y2 (xor) o Y2 (xiT) , 

is actually a fourth order symplectic integrator, that moreover is symmetric, i.e. it 
has exact time reversibility. 

One can iterate this construction and find new explicit coefficients, yi = 2-2'^/^ 
and yo — —2^/^yi, such that the composition 

(21) l6(r) = Yi (yiT) o ^ (j/qt) o Y4 (yir) , 

provides a sixth order symmetric symplectic integrator. This idea can be iterated to 
construct symmetric symplectic integrators with arbitrarily high order. The main 
drawback is that the number of involved terms increases very fast, thus one has to 
choose a suitable compromise between the required precision in term of preserved 
quantities, i.e. energy and mass, and the CPU time available. 

In the next sections we will show that 1^4 and Yq exhibit very good energy and 
mass preservation properties even for relatively large integration time steps, they are 
composed by relatively few terms and moreover because they are explicit methods. 
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they are relatively fast. They provide thus very powerful and fast methods to 
numerically analyze the DNLS in the asymptotic limit of large time spans. 

Remark 3.3. Let us observe that our method straightforwardly applies to general- 
ized DNLS [3] 



HgDNLs{p,q) = 



^ ^ {pi+i - pi)im+i - qi) f •,a+2 
(piqi) 



1=1 



being a a positive parameter. In fact setting Aa{p,q) ~ ^^iLiiPili)"^^^ '^''^d, ob- 
serving that Ci = piqi is still a first integral for the flow of A^, we directly obtain 
for exp(TiA„) 



V/e {!,..., iV} 



While the flow associated to the B part remains unchanged. 

More generally our method can be directly applied whenever we replace the A-part 
of the Hamiltonian function with a new one, A'(p,q), for which the map exp(rL^') 
can be computed explicitely. 

3.4. Preserved quantities. By construction the symplectic scheme (t) will 
preserve the energy of the systems with an error O (t^™) , that is independent of the 
relative weight of the functions A and B used to decompose the Hamilton function, 
this is the reason why this method is more suitable than the one proposed in [S] 
where the error is also a function of the relative weights. 

The aim of this section, is to prove that the symplectic schemes Y2m preserve 
other relevant quantities of the DNLS system ([3]). 



3.4.1. Preservation of the first integral I{p, q) = piqi . To prove that our method 
preserves the first integral / let us start by consider the action of the map exp (tLb) 
on the function / — ^iPiqi. Starting by its very first definition (I13p we get: 

d i ^ 

^/|flowi3 = -j^^[{pi+i-2pi+pi^i)qi-pi{qi+i-2qi+qi^i)] 
1=1 

. N 

= -~j^^ipi+iii +PI-111 -piqi+1 -pm-i) = 0, 
1=1 

where the last equality follows by using the boundary conditions. Thus the flow of 
B preserves /. 

On the other hand by the definition of the map exp {La) given by (jlip one 
straightforwardly get: 

Pl{t)qi{t) = pi{Q)qi{Q) €{!,..., TV}, 

and thus also the flow induced by A preserves /. 

We can thus conclude that the composition of the maps exp {La) and exp {Lb) 
preserves the first integral / and so does any symplectic scheme Y2m- 
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3.4.2. Preservation of the relation p{t) — q{t). Let us define the complex vector 
A{t) = p{t) — q{t). Using tlie fact ttiat M is real, the time evolution of A under 
the action of B is given by ([T5|) : 



from which one gets: 



A(t) = e--^*"A(0). 



Thus if by assumption A(0) = then we get A(t) — for all t. 

For the flow of A we use once again the explicit map (fill) and the fact that 
^ = (f(0) implies that C/ = pi{0)qi{Q) = |g/(0)p is real for alH e {1, . . . , A^}, thus 

^ = e^^^'*^ = e'^^'*gi(0) = qi{t) . 

Once again, we can thus conclude that the composition of the maps exp {La) 
and exp [Lb) preserves the first integral / and so it does any symplectic scheme 

3.4.3. Other preserved quantities. The symplectic integrators we proposed preserve 
other quantities such as the norm of the canonical variables, namely and |(?|^, 
defined by the complex scalar product ^< p,p >. 

Let us first observe that this statement holds for the flow of B. For instance in 
the case of the p variable one has: 

d ^ ^ ^ ^ i _ i ^ 

— <p,p> = < p,p> + <p,p>=Eq. ^= < Mp,p> <p,Mp> 

= ^ < Mp,p> < M''p,p>=0, 

where in the last step we used the facts that M is real, thus = M'^ , and 
moreover M'^ = M. 

It remains to check the behavior under A. But using the definition we get: 

I I I 

Along a very similar way we can prove the result for q. We can thus conclude 
that the symplectic integrators preserve the norm of the complex vectors p and 



4. Results 

The aim of this section is to present the numerical results obtained using our 
high order symplectic schemes. We fixed initial conditions as in |101 [S] to have 
a good testbed to compare our results with other ones available in the literature, 
hence we define: 

(22) 9,(0) =pKO) = a(l-ecos(teO) , 

where xi = —L/2 + {I — l)/i, h = L/N , I G {1, . . . , N}, namely we are considering 
perturbation of a spatially uniform plane wave invariant under the phase flows. We 
used several values for N while the remaining parameters have been fixed jlOl [5] to 
e = 10-2^ b = 2ti/L and L = 2^/2ti. 
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Let US stress that because of the high accuracy of most of the presented results, 
mainly in term of energy preservation, we performed our numerical simulation using 
quadruple precision Fortran. 

The first result reported in Fig. [T] shows the preservation of the energy for one 
given orbit with the above initial conditions using the symplectic schemes 
for m e {1,2,3,4}. We can observe that the relative energy loss, AE{t) = \E{t) — 
i?(0)|/|i?(0)|, fluctuates in times but doesn't grow on a quite large time span [0, 10"'] 
even using a relatively large time step r = 10^^. Moreover we can remark that 
already with Y4 the relative energy loss, is well below 10"^^°. On the other hand Yg 
allows to reach values of the order of 10"^". 

The preservation of the other quantities is even better; concerning the mass, 
I{p,q), we can find that the relative error, A/(i) = — /(0)|/|/(0)|, assumes 
values well below 10^^^ for all integrator schemes we used Y2m (data not reported); 
while the relation p{t) — q(t) is preserved up to the machine precision, namely 
~ XQ-34 ^(ia,ta not reported). 

In the inset of Fig. [T] we report the comparison with the non-symplectic integrator 
Runge-Kutta 4 using r — 10^^; we can clearly see the inefficiency of the latter 
method even using a time step smaller than the one used for the symplectic schemes, 
in fact the relative energy loss becomes larger than 10^ already at i ~ 350, so no 
longer comparison are possible. Using Runge-Kutta 4 the relative energy loss 
grows in times, the slower is the time step, but still growing, hence one can reach 
the precisions obtained by a symplectic integrator only using very small time steps 
T or integrating over relatively short time spans. For instance one can achieve a 
relative energy loss of ~ 10~^ using RK4{10~'^) only on a time span [0,^ 350], 
while using l2(10~^) we can get the same error but on [0, 10'']. On the time span 
[0, ^ 350] and using the time step r = 10^'', Runge-Kutta 4 achieves a relative loss 
for the mass of the order of 10"^ while the conjugacy relation p(t) = q{t) is not 
preserved any more. 

This is the main reason of the poor properties obtained using Runge-Kutta 4 
with r — 0.01; in fact if we modify the integration scheme by forcing the vectors 
p{t) and q{t) to satisfy the conjugation relation at each time step, we obtain a 
method, hereby called modified Runge-Kutta 4, RK^"''-, that exhibits improved 
energy preservation properties (see Fig. [1]). Let us observe that this is method is 
not symplectic, as one can clearly conclude from the increasing trend in the relative 
energy variation presented in the Figure. The method allows to reach large time 
spans, [0, 10''], but the goodness of the method, measured in terms of relative energy 
variation, is worse than 14 and slightly better than Y2. Let us finally observe that 
on the same large time span and still using r = 0.01, the mass of the system is 
conserved up to a factor 1.5 10~^. 

The next study will concern the dependence of the maximum relative energy 
loss as a function of the time step, by integrating several solutions using Y2m (t), 
m S {1,2,3,4,5}, over a large time span [0,10'']; let us observe, once again, that 
because the relative energy fluctuates around the flxed initial value, the same results 
hold for arbitrarily longer time spans. Result reported in Fig.[2]shows the computed 
numerical accuracy of Y2m (t) as a function of the time step r. Let us observe 
that the use of Yio, respectively of Yg, with integration steps r smaller than ~ 
2. 10~^, respectively ^ 8. 10"'', will produce a maximum relative energy loss below 
the quadruple machine precision, that is why we limited our simulations to these 
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100 200 300 400 



2000 4000 6000 8000 10000 

time 



Figure 1 . Time evolution of relative energy. Semilog plot of the 
relative energy loss, AE{t) = \E{t) - £'(0)|/|£'(0)|, as a function 
of time, for one orbit with initial conditions (|22t . = 4 and a 
time step r = 0.01. Using, from top to bottom, the integrators 
Y2 (blue), RK2'°'^ (magenta), Y4 (red), (black) and Fg (green). 
Inset : using RK4 with time step r = 0.001. 

values. Straight lines reported in Fig. [5] represent linear best fits logj^Q max | AE'j = 
a logj^o + /5j whose coefficients a and /3 are reported in Table [I] and numerically 
confirm the accuracy of the integrators Y2m- 

Let us remark that a similar result holds for larger values of N. The main 
difference being that in this case large time steps are prevented to be used because 
of a decrease in the performances of all integrator schemes, F2m('''), mainly because 
of the energy preservation. This fact has been already observed [5] and relies on a 
stability issue of the integrators, symplectic and non-symplectic ones, that imposes 
a constrains tN « c, for some positive constant c. 

Our last remark concerns the speed of the numerical methods (t) ■ In Fig. [3] 
we report the CPU time used to integrate orbits with initial conditions (1^^ and 

= 4 using l2„i(T), as a function of the time step r. For Y^ and Yiq we limited 
the analysis to time steps whose maximum relative energy loss is larger than the 
machine precision using quadruple precision Fortran (see Fig. [5] and discussion 
therein). From these data we clearly conclude that the CPU time increases as l/r 
for a fixed integrator scheme l2m(T); on the other hand, for fixed r, the CPU time 
increases as a exponential of the integrator order, roughly as 3™. This is because, 
as already mentioned, the Yoshida scheme does not get the optimal number of 
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10~^ 10^ 10^ 10 



Figure 2. Accuracy of the integrators as a function of t. We 
report maxtgjg ]^o4] |Ai?(i)| for orbits with initial conditions given 
by (j22p and = 4, numerically integrated over the time span 
of [0,10*] with O Y2 (black), □ Y4 (blue), Yq (green), A Fg 
(red), V Yio (cyan), empty ★ RK^°'^ (grey) and ★ RKi (grey) 
but over a small time span [0,200]. Straight lines are linear best 
fits logj^o max I Ai?| = alogj^gr + /3, the coefficients a and /3 are 
reported in Table [TJ 



integrator 


a 


/3 


Y2 


1.997 ±0.012 


-1.42 ±0.03 


I4 


3.996 ±0.009 


-2.271 ±0.023 


Ye 


5.97 ±0.06 


-2.23 ±0.15 


Ys 


7.97 ±0.07 


-1.86 ±0.16 


Yw 


9.88 ±0.17 


-1.7±0.3 


RKi 


3.997 ±0.001 


-0.716 ±0.006 




4.40 ±0.11 


0.82 ±0.29 



Table 1. Numerical coefficients of the linear best fits 
log max \ AE\ — a logj^Q r + /? reported in Fig. [21 



products exp{LA) and exp{LB)- On one hand Yoshida already proposed [12] an 
improved version to tackle this difficulty, whose results are symplectic schemes with 
fewer compositions ([7]). Because the computations to obtain the coefficients {cj,dj) 
become rapidly cumbersome we limit ourselves to study the cases of sixth order. 
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Yg"''*, and eighth order, Yg"^*. Results reported m Fig. [3] show that Ye{T) needs 
almost 1.5 times more CPU time than Yq''\t), while Ys{t) about 2.2 times more 
CPU time than Ys^'^t). 

On the other hand in practical applications one should choose the time step r 
and the integration order 2m that produce the best balance between the required 
precision, say maximum relative energy loss, and the CPU time needed. For in- 
stance from Fig. [2] we clearly see that using Y4 and r ~ 10"'^, we can ensure a 
maximum relative energy loss of the order of 10^^^, to get the same precision using 
Ye one can use a time step ten times larger, t ^ 10^^ and thus (see Fig. [S]) the re- 
quired CPU time will be smaller with 16(10"^) than with l4(10~"^). Requiring the 
same precision, using instead Ig, will need a time step "only" twenty times larger 
and thus the CPU time wiU increase: Ys{2. 10"^) > Yq^IQ-'^). 

The Runge-Kutta 4 requires a CPU time larger than Y4 using the same time 
step, roughly of the same order of Yg"^*. On the other hand the modified Runge- 
Kutta scheme is faster than Y4 because it has to compute only half of the vector 
field. 

The dependence of the CPU time on the discretization parameter h, hence on 
N, is more crucial once we need to use very large N and/or a very large number of 
orbits. Because in the present work we were not interested in the optimality of the 
method, we computed "naively" the map exp(LB), namely using a vector-matrix 
product whose cost is O (-/V^) . On the other hand we can easily improve this part by 
considering the strong similarity with the map exp(Ls) and the Fourier transform 
of the vectors {p,q) (see Proposition 13. II and Eq. (HH])) and thus use, instead of the 
vector-matrix product, a Fast Fourier-like method to speed up the computations. 



5. Conclusions 

In this paper we presented a family of high order, explicit, symplectic integra- 
tions schemes adapted to the study of the DNLS. Despite DNLS has been studied 
numerically since long time, this is the first time that such a high precision can be 
achieved using relatively large time steps. Besides the very good energy preserva- 
tion properties of the above introduced methods, we also obtained an almost exact 
preservation of the other first integral, the mass of the system, and of the conjugacy 
relation. Because the integrators we constructed are explicit ones, they result very 
fast. 

For all these reasons we believe that such accurate numerical schemes could be 
very useful to test several physical hypotheses concerning the asymptotic regimes 
of the DNLS, for instance the existence and stability of breathers and the regimes 
with negative temperature. 
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Figure 3. CPU Time needed by Y2,n{T) as a function ofr. For 
a fixed time span of [0, 10'*] we report the CPU time needed by 

(t) to integrate orbits with, initial conditions given by 
and = 4. Symbols are : Q Y2 (black), □ (blue), Yq 
(green), empty Yq^* (green), A Yg (red), empty A Fg°^* (red), 
V ^10 (cyan), empty ★ RK]^°'^ (magenta) and ★ RK4 using (ma- 
genta) . 
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